Momentum-dependent pseudogaps in the half-filled two-dimensional Hubbard model 
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We compute unbiased spectral functions of the two-dimensional Hubbard model by extrapolating 
Green functions, obtained from determinantal quantum Monte Carlo simulations, to the thermody- 
namic and continuous time limits. Our results clearly resolve the pseudogap at weak to intermediate 
coupling, originating from a momentum selective opening of the charge gap. A characteristic pseudo- 
gap temperature T* , determined consistently from the spectra and from the momentum dependence 
of the imaginary-time Green functions, is found to match the dynamical mean-field critical temper- 
ature, below which antiferromagnetic fluctuations become dominant. Our results identify a regime 
where pseudogap physics is within reach of experiments with cold fermions on optical lattices. 

PACS numbers: 71.10.Fd, 71.27. +a, 71.30.+h, 74.72.-h 



I. INTRODUCTION 

A peculiar feature of (undcrdoped) high-T c supercon- 
ductors is the coupling of antiferromagnetic fluctuations 
to charge degrees of freedom, which leads to a strong mo- 
mentum dependence of the spectral functions. In partic- 
ular, it induces pseudogaps in the normal state, i.e., a 
suppression of the density of states at the Fermi energy, 
which can be probed using (angular resolved) photo- 
emission, inverse photoemission, and related techniques. 
The pseudogap shares the d- wave symmetry with the or- 
der parameter of the superconducting phases occurring 
at low temperatures and near optimal doping 

Pseudogap physics can also be expected in the undoped 
Hubbard model. In the absence of electronic correlations, 
the tight binding model is characterized by a coherence 
temperature T co h, set by the bandwidth W. For weak 
Hubbard interaction U < W, antiferromagnetic spin fluc- 
tuations, with energy scale T sp i n , will develop below the 
coherence temperature. Hence, the temperature window 
Tjv < T < T S pi„ < T co h is characterized by a metallic 
state coupled to antiferromagnetic fluctuations, which 
sets the stage for pseudogap physics. Here Tjv is the 
Neel temperature, at or below which long range order 
generates a full gap in the presence of perfect nesting (in 
dimensions d > 2, with Tjv = in d = 2). 

Theorists have tried to verify this scenario on the 
basis of numerical simulations and to compute reliable 
spectra for decades. Direct simulations can only be 
performed for clusters of finite extent, usually employ- 
ing periodic boundary conditions. Early determinantal 
quantum Monte Carlo^ (DQMC) studies at moderately 
weak coupling (U/t — 4) led to spectra with significant 
low-temperature pseudogap features only for small clus- 
ter sizes; thus, pseudogaps in the undoped 2d Hubbard 
model were regarded as pure finite-size (FS) artifacts^ 
Later studies at similar coupling strengths^— found 
pseudogaps also at large cluster sizes, but did not al- 
low for quantitative predictions in the thermodynamic 
limit. A recent study using the dynamical vertex approx- 
imation (DrA) observed reentrant behavior incompatible 



with the earlier resultsAi 

A central limitation of all previous numerical pseudo- 
gap studies is that results for different cluster sizes (e.g. 
in DQMC simulations) were compared only at fixed tem- 
peratures and at the level of spectral functions. With in- 
creasing cluster size, these show diverse effects: shifts of 
spectral peaks, transfer of spectral weight, and the open- 
ing or closing of gaps. A direct pointwise extrapolation 
of these positive semidefinite and normalized functions 
is clearly impossible. In fact, we are not aware of any 
published attempts of deriving spectral properties in the 
thermodynamic limit from DQMC data in any context. 

In this paper, we present (i) the local spectral func- 
tion, (ii) momentum-resolved spectral functions at high- 
symmetry points, and (iii) momentum-resolved spectral 
functions along high-symmetry lines of the Brillouin zone 
in the thermodynamic limit. All results are based on 
FS extrapolations of imaginary-time Green functions, ob- 
tained from DQMC, with subsequent analytic continua- 
tion to the real axis using the maximum entropy method 
(MEM)^ and, in case (iii), a Fourier fit of the momen- 
tum dependence. The final results are free of significant 
systematic errors and represent the thermodynamic and 
continuous time limits in an unbiased way. 

Thereby, we can not only unambiguously confirm the 
pseudogap scenario and study the nodal-antinodal di- 
chotomy in unprecedented detail, but also explore the 
temperature dependence of the pseudogap opening and 
identify a characteristic temperature T*. At weak to in- 
termediate couplings, T* tracks the onset of short-ranged 
magnetic fluctuations, and is equally shown to compare 
remarkably well with the dynamical mean-field critical 
temperature for antiferromagnetic long-range order. 

In Sec. [Ill we introduce the model, set up our notation, 
characterize the established methods (DQMC, MEM) 
underlying our approach, and specify our implementa- 
tions. The new methods for eliminating systematic biases 
from Green function and spectra are, then, presented in 
Sec. IIII1 first for the DQMC Trotter error, then for finite- 
size effects. Our main results are discussed in Sec. IIV| 
starting with pseudogap features in the spectral func- 
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tions for the "nodal" and "antinodal" — high-symmetry k 
points on the Fermi surface and their evolution as a func- 
tion of temperature. We then show, with continuous mo- 
mentum resolution, how the pseudogap evolves through- 
out the Brillouin zone (BZ) and discuss non-Fermi liquid 
physics that is not accessible by conventional methods. 
Finally, we determine the characteristic pseudogap tem- 
perature T* for U < W and relate it to spin correlation 
functions and other characteristic temperature scales of 
the model. 



II. MODEL AND CONVENTIONAL METHODS 

A. Hubbard model 

Our starting point is the single-band Hubbard Hamil- 
tonian with nearest-neighbor hopping t on a square lat- 
tice (with unit lattice spacing a = 1): 

H = Ho + U^n^hn (!) 

i 

Ho = -tj^ Csv = E £kftk - ( 2 ) 

(ij),<r k,o- 
£ k = -2t[cos(k x ) + cos(ky)} (3) 

Here, c- (cj ) are annihilation (creation) operators for a 
fermion with spin a G {1\ i} a t site i; hi a = c\ a c ia . In 
the following, the energy scale will be set by t = 1. 

At half filling n = (n^ + n^) = 1 and in the noninter- 
acting limit {7 = 0, the occupied momentum states form 
a square (dark shaded) within the square Brillouin zone 
illustrated in Fig. [TJa) (in the thermodynamic limit), 
which implies a perfect nesting instability: The Fermi 
surface (gray line) transforms into itself when shifted by 
the antiferromagnetic wave vector (tt,tt). As a conse- 
quence, any finite interaction U > drives this model to 
long-range antiferromagnetic order (only) in the ground 
state. 

While a conventional notation is well established for 
the center T and the corner M of the BZ of the square 
lattice as well as for the antinodal X point, this seems 
not to be the case for the nodal point; in Fig. Q] and in 
the following, we denote this midpoint of TM as M'. 

B. Determinantal quantum Monte Carlo algorithm 

The Hubbard model, Eq. (1), is solved in this work 
for clusters with a finite number N of sites [implying 
a discrete momentum grid, see Fig. [TJb)] at finite tem- 
peratures T using the DQMC algorithm developed by 
Blankenbecler, Scalapino, and Sugar j£ with modifica- 
tions by Hirsch.— It is based on (i) a uniform discretiza- 
tion of the imaginary-time interval < r < j3 [with 
/3 = l/(fcsT)], occurring in the path integral, into A 
time slices of width At = (3/ A, (ii) a Trotter decoupling 




FIG. 1. (Color online) Brillouin zone (BZ) of the square lat- 
tice: (a) full BZ with Fermi surface (gray line) and occupied 
momenta in thermodynamic limit (dark shaded) at U — 0; 
letters denote high-symmetry points; color lines and arrows 
indicate the path used in Fig. [9] (b) Irreducible wedge of 
BZ with momenta occurring in finite-size clusters of linear 
dimension L = 8, 10, 12 with periodic boundary conditions. 

of kinetic and interaction terms, and (iii) a Hubbard- 
Stratonovich transformation which replaces the electron- 
electron interaction at each time slice and site by the 
coupling of the electrons to a binary auxiliary field. Ex- 
pectation values are obtained by Monte Carlo importance 
sampling of field configurations, with weights given by 
a product of two determinants for the two spin compo- 
nents. In the particle-hole symmetric case considered in 
this study, this product is always positive, i.e., the sign 
problem is absent. The numerical effort scales as f3N 3 . 
A detailed review of the algorithm used can be found in 
Ref . EE 

In this work, we obtained imaginary-time Green func- 
tions and spin correlation functions between each pair of 
sites by applying the DQMC method to square lattice 
clusters L x L of the linear size L = 8, 10, 12, 14, 16 with 
periodic boundary conditions, using a set of Trotter dis- 
cretizations with 0.1 < At < 0.42 and typically 50 bins 
with 5000 sweeps over the auxiliary field each. For the 
largest systems, individual runs took about a month of 
computer time; up to five such runs were averaged over 
in order to reduce error bars. This resulted in typical 
statistical errors in the (finite-size) Green functions of 
O(10 -4 ). Note that the DQMC scaling with L e makes 
it difficult to access much larger system sizes directly: 
L = 20 (L = 32) would increase the effort by a factor of 
about 4 (64) compared to L = 16. Local properties were 
averaged over all sites, momentum (k) dependent proper- 
ties were obtained by Fourier transforms of the real-space 
measurements. 



C. Maximum entropy method 

Since DQMC calculations provide Green functions G 
(and correlation functions) only at imaginary times, their 
interpretation as dynamical information requires an an- 
alytic continuation to the real axis. Specifically, one has 



to invert relations of the form 

/oo e~ TU 

where A(ui) = —ImG(uj)/ir is the corresponding spectral 
function. This is an ill-posed problem, as the exponen- 
tial kernel suppresses the impact of features in A(u>) at 
large \u\ on G(t); in the DQMC context, further com- 
plications arise from the fact that G is only measured 
on the discrete imaginary-time grid {r; = ZAt}^ 1 . The 
MEM finds the most probable spectrum, given the data 
Gi ± AG/, by balancing the misfit of a given candidate 
spectrum (measured by the corresponding y 2 ) with an 
entropy constraint which favors smooth spectra i 12 ' 16 In 
our implementation, the resulting minimization prob- 
lem is solved deterministically using a Newton scheme 
in the singular space of the kernel. Its application both 
to DQMC raw data for local and k dependent Green 
functions and to Green functions obtained from Trotter 
and/or FS extrapolations always resulted in reliable and 
consistent maximum entropy spectra. 



III. EXTRACTION OF UNBIASED SPECTRA 
IN THE THERMODYNAMIC LIMIT 

What sets our main results, to be presented in Scc. lIVl 
apart from earlier work, is their direct relevance in the 
thermodynamic limit, i.e., the absence of significant bias. 
We now specify our methodology for eliminating Trotter 
and finite-size errors from Green functions and establish 
its accuracy and reliability on the level of Green functions 
and spectra. 



A. Trotter errors and extrapolation At — > 

As discussed in subsection III B[ the DQMC method 
decouples electronic interactions (and evaluates, e.g., 
Green functions) at the cost of introducing an artificial 
imaginary-time discretization At, which implies an un- 
physical bias in all DQMC estimates of obscrvablcs. In 
the absence of phase transitions, DQMC raw results are 
expected (and observed) to depend smoothly on At, in 
the form of a power series; for some static observables, 
such as the total energy, it is easy to proved polynomial 
dependence on At 2 . 

The effects of the Trotter discretization on the 
imaginary-time Green function G(t) are illustrated in 
Fig.[2ja): (i) each of the raw data sets (symbols) lives on 
a different r grid; (ii) at fixed values of r, the data (or 
a linear interpolation - dashed/dotted lines) is shifted to 
smaller absolute values at larger At. Obviously, unbi- 
ased results (for a fixed cluster size L in real space) can 
only be expected after an extrapolation of At — > 0. On 
the other hand, such an extrapolation is not possible lo- 
cally, i.e., at fixed t, but requires a global approach that 
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FIG. 2. (Color online) Impact of Trotter discretization and 
extrapolation of At — > 0: (a) for local imaginary-time Green 
functions G(t) at U = 4, T = 0.20, L = 8 in the range 
< t < fi/2 — 2.5 (inset: magnified view for r > 1.7); 
(b) corresponding local spectral functions A(u). (c) for the 
difference between nodal and antinodal Green functions (cf. 
subsection IIV C[) versus temperature. 

can use input from DQMC raw data at all discretizations 
At for each imaginary time t of interest. 

The black solid lines in Fig. EJJa) represent the result 
of a multigrid procedure, originally developed in the con- 
text of the Hirsch-Fye quantum Monte Carlo method for 
the Anderson impurity model.— The multigrid method 
is based on the fact that "reference" Green functions 
G re f(T) with sufficiently accurate asymptotics at t — > 
(and t — > j3), in particular for the curvature, can eas- 
ily be derived from weak-coupling expansions (or, alter- 
natively, from the "best" QMC data via MEM); con- 
sequently, the difference between the measured Green 
functions Gat(t"i) and the reference G re f(Ti) can be ad- 
equately represented by a natural cubic spline for each 
value of At; all of these splines can, then, be evaluated 
on a common (fine) grid. For this transformed data, we 
find a nearly linear dependence on At 2 (plus a small cur- 
vature) , so that pointwise extrapolations At — > are re- 
liable and accurate. At the level of G, the shift of the un- 
biased result [black solid line in Fig. H{a)] of about 10 -3 
compared to the best raw data [at At = 0.1, squares and 
dash-dotted line in Fig. EJa)] is still significant. 

This is no longer true on the level of spectra, shown in 
Fig.[2jb), due to the intrinsic complications of MEM: The 
results for At = 0.1 agree within accuracy with the unbi- 
ased spectrum. Thus, we may conclude that At = 0.1 is 
"good enough" for spectral data (at U = 4) and that an 
elimination of the Trotter error is not necessary for reduc- 
ing unphysical bias below significance. At the same time, 
the smooth consistent evolution of the spectra with At 
confirms our MEM procedure both for the DQMC raw 
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FIG. 3. (Color online) Imaginary-time Green functions (for 
selected values of r) versus squared inverse linear dimension 
(empty symbols) and least-squares extrapolation L — > oo 
(lines and large full symbols) at high-symmetry momentum 
points: (a) at antinodal X point, (b) at nodal M' point (small 
symbols for L = 10, 14: from Fourier fit as shown in Fig. 0). 

data and for extrapolated Green functions. 

Even smaller Trotter errors than observed in Fig. [DJa) 
can be expected for differences of Green functions, due 
to error cancellation. Indeed, the scalar pscudogap mea- 
sure |Gm' — Gx\, to be introduced in subsection IIV CI 
is impacted significantly by Trotter errors only for large 
discretizations; the bias become negligible for At < 0.1, 
as shown in Fig. EJc). Therefore an explicit elimination 
of this error is, again, not necessary. 



B. Finite-size extrapolations of local spectra and 
on high-symmetry k points in the BZ 

A FS extrapolation of local properties or k resolved 
properties at high-symmetry points is relatively straight- 
forward: One accumulates raw data at various values of 
the linear extent L and then extrapolates using polyno- 
mial least-square fits in 1/L 2 . In the case of imaginary- 
time Green functions, independent extrapolations have 
to be performed for each value of r (on the grid with 
spacing At in the case of DQMC raw data or the grid 
chosen in the extrapolation At — > discussed in the pre- 
vious section). As shown in Fig.[3]for U = 4, T = 0.2, the 
Green function depends on system size quite significantly 
at generic imaginary times [except for the limits r — > 
or, equivalently, t — 5- (3 (not shown)] both at the antin- 
odal (a) and nodal (b) momentum points. At the same 
time, the dependence is quite regular so that least-square 
extrapolations (lines) can be restricted to low orders. 

Obviously, this "local" procedure can only include lat- 
tice sizes for which the k point under consideration exists 
[cf. Fig.QJb)]; for the antinodal point, this requirement is 
fulfilled for all even values of L, while the nodal M' point 
is only present if L is a multiple of 4 (which restricts the 
set to L = 8, 12, 16 in our study). Still, as seen in Fig. 
[3jb), the extrapolation is reliable even with only three 
data points (per fit), as the dependence is almost per- 
fectly linear At the same time, the FS extrapolation is 
particularly important at the nodal M' point, as FS ef- 
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FIG. 4. (Color online) Upper row: imaginary-time DQMC 
Green functions (U = 4, T = 0.2, At = 0.1) at antinodal 
(a) and nodal (b) points, respectively, for finite-size clusters 
(open symbols and broken lines) plus extrapolated (cf. Fig. [3]) 
results in thermodynamic limit (filled circles and solid lines). 
Lower row: corresponding spectral functions. 



fects are much stronger than in the antinodal case [shown 
in Fig. Ha)]. 

Note that 4x4 clusters (with periodic boundary con- 
ditions) have a special symmetry: They have the same 
topology asa2x2x2x2 hypercube with open boundary 
conditions; as a consequence, the next-nearest neighbors 
along one of the axes and the ones along the diagonal be- 
come equivalent, which implies that the X and M' points 
are identical in momentum space at L = 4. In order to 
avoid the associated extra bias we exclude this system 
size and consider only lattices with L > 8 in this study. 

The full resulting Green functions in the thermody- 
namic limit are shown as solid lines in Fig. U for the 
antinodal (a) and nodal (b) k points, respectively, to- 
gether with their finite-size equivalents (dashed and dot- 
ted lines). We see, again, that FS effects are much more 
prominent at M' [note the different scales in the insets of 
Fig. BJa) andlUJb)]. The effect is even much stronger on 
the level of the corresponding spectra, shown in Fig.[4jc) 
and0Jd), respectively: In an 8 x 8 system (dashed-dotted 
line), nodal and antinodal spectra are qualitatively very 
similar, with a clear pscudogap feature, and differ mainly 
in peak height (at \u\ 0.7); at k = X, the spectrum 
remains nearly unchanged at larger system sizes and in 
the thermodynamic limit. At k = M', in contrast, the 
pscudogap dip shrinks significantly for larger systems and 
is completely lost in the thermodynamic limit, where a 
quasiparticle shape appears. This shows that essential 
pseudogap physics, with a nodal-antinodal dichotomy, is 
really a property of the thermodynamic limit and that 
the bias inherent in finite-size systems dangerously dis- 
torts the physical picture. 
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FIG. 5. (Color online) Example of finite-size extrapolation of 
imaginary-time Green functions along high-symmetry lines in 
the BZ, here for the path X -»■ M' [cf. Fig. Ufa) ] and r = /3/2: 

(a) difference Green functions with respect to Gx (symbols), 
fitted with a Fourier series (broken lines) and final result of the 
extrapolation to the thermodynamic limit (black solid line), 

(b) extrapolation of the corresponding Fourier coefficients. 
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C. FS extrapolations of spectra along 
high-symmetry momenta in the BZ 

The elimination of FS errors at generic momenta re- 
quires "global" extrapolations that involve some kind of 
functional fitting procedures in momentum space. For 
momenta along high-symmetry lines, these fits have the 
form of Fourier series which may be adapted in order 
to take all symmetries into account. In the following, we 
will illustrate the algorithm for the most important path, 
the irreducible portion XM' of the noninteracting Fermi 
surface. This path can be parametrized as 

k x = (2 - k) tt/2; fc y = K7r/2; « € [0, 1] ; 

then k = corresponds to X e (tt, 0) while k = 1 corre- 
sponds to M' ee (tt/2, tt/2). At particle-hole symmetry, 
all functions / have to be symmetric with respect to both 
end points, which implies that they can be represented 
in the form 

oo 

f(n) = ao + a n sin 2 (nK7r/2) 

n=l 

with coefficients a n . We have chosen to fit the difference 
of the Green function for each k (along the line) with re- 
spect to the antinodal Green function (corresponding to 
k = 0); this implies that the zeroth-order coefficient van- 
ishes exactly. The symbols in Fig. [3a) represent DQMC 
data for the difference Green functions at r = (3/2 = 2.5; 
evidently their interpolation using the above Fourier se- 
ries up to third order (dashed/dotted lines) works quite 
well. Furthermore, the associated Fourier coefficients de- 
pend very regularly (i.e., almost perfectly linearly) on 
1/i 2 , as seen in Fig. E^b), and decay exponentially as 
a function of order. Consequently, an extrapolation to 
the thermodynamic limit is possible on the level of the 
coefficients (using a least-squares fit) with high preci- 
sion; the extrapolated coefficients yield a reliable esti- 
mate Gk(T = 2.5) for all k along the path [solid line in 



FIG. 6. (Color online) Evolution of the DQMC spectral 
functions with temperature at U = 4 for finite clusters with 
8 < L < 16 and in the thermodynamic limit: (a) at the 
antinodal point [k = X = (tt, 0)], (b) at the nodal point 
[k = M' = (tt/2, tt/2)]. 

Fig- E^a)]. This procedure has to be performed indepen- 
dently for each value of r; spectra can then be obtained 
using MEM on an arbitrarily dense k grid. Similar proce- 
dures were employed separately for each high-symmetry 
line indicated in Fig. []Ja). 



IV. RESULTS 

A. Pseudogap signatures at nodal and antinodal k 
points 

Let us, first, turn to the antinodal and nodal spectra 
shown in Fig. [SJa) and [S{b), respectively. At the ele- 
vated temperature T = 0.50 (dotted lines) the spectra 
have quasiparticlc (QP) shape at all system sizes and for 
both momentum points. FS effects are negligible: Even 
the spectra of the smallest systems considered (8x8, left 
column) do not deviate visibly from those in the ther- 
modynamic limit (right column); also the momentum de- 
pendence along the Fermi surface is minimal at T = 0.50, 
with about 20% larger peak height at the nodal M' point. 

In the 8x8 case (left column), the largest system size 
fully considered in previous studies, a pseudogap dip ap- 
pears almost simultaneously at T = 0.28 and T = 0.24 
(dashed lines) at the X and M' points, respectively, and 
quickly deepens to an almost complete gap at T = 0.18 
(dashed-dotted line). Given only this data, one would 
conclude that any momentum dependences beyond the 
free dispersion are inessential, i.e., that the physics might 
be in reach of theories with a momentum-independent 
self-energy [in particular, the dynamical mean-field the- 
ory (DMFT)]. However, this picture is distorted by finite- 
size effects and far from the truth: In the thermodynamic 



FIG. 7. (Color online) Scalar measure of pseudogap strength 
(see text) versus temperature: (a) at antinodal, (b) at nodal 
point. The nodal-antinodal dichotomy is fully apparent only 
in the thermodynamic limit (solid lines). 

limit (right column in Fig. [5]) , the antinodal spectra have 
QP shape only for T > 0.28; at T = 0.24, a slight dip 
emerges at ui = which develops to a significant PG at 
T = 0.20 and an almost complete gap at T = 0.183 
In contrast, the nodal spectrum retains QP form down 
to T = 0.20 (while even the 16 x 16 system shows a 
PG dip at this temperature), before a PG emerges at 
T = 0.18. Thus, the FS extrapolation detailed above 
is really essential for fully resolving the nodal-antinodal 
dichotomy, which is at the heart of PG physics at fi- 
nite temperatures. Only in the limit T — > 0, i.e., in the 
presence of long-range order, one expects a DMFT-likc 
picture to become valid (again, as for high T) with fully 
gapped spectra all along the Fermi surface. 

This implies that finite-size effects should mainly have 
two consequences on the level of spectra: (i) shift charac- 
teristic PG temperatures upwards, (ii) dilute the nodal- 
antinodal dichotomy in the vicinity of these characteristic 
temperatures. 

Characteristic PG temperature T* - As the opening 
of the PG is not a thermodynamic phase transition, it 
lacks a unique critical temperature. It is still useful (and 
common)**"* to define a characteristic PG temperature 
T* , for comparison with other temperature or energy 
scales of the system. An obvious choice of the required 
scalar PG measure is a dip in the spectral function. We 
specify this "pseudogap strength" by the reduction of 
spectral weight at u = 0, compared to the maximum: 

rk = 1 — Ak_(u> = 0)/ maxAk(w) , 

as shown for the (anti)nodal points in Fig. [3 This rep- 
resentation reveals that the onset of the PG is slow only 
at k = X: As soon as rx ~ 0.5, t°m' jumps to the full 
value within a narrow temperature range AT w 0.03. 
The results in the thermodynamic limit (filled circles) 
can be fitted with a Fermi function form (solid lines); us- 
ing their inflection points yields T x « 0.20, Tfc, « 0.18. 
Note that, again, the FS effects are much stronger at 
k = M' than at k = X. 

Comparison with the literature - In a pioneering study, 
Huscroft et a/3 had obtained first bounds on the FS 
errors in DQMC spectra by complementing DQMC re- 
sults for N < 64 sites with those of the dynamical clus- 
ter approximation (DCA) employing N k patches in the 
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FIG. 8. (Color online) Spectral functions of the half-filled 
Hubbard model at the antinodal point [k = X = (it, 0)] and 
for weak coupling U = 5.2: (a) Unbiased spectrum at T — 
0.20 (solid line), in comparison with earlier DCA and finite- 
size DQMC results.— (b) DQMC spectra for 12 x 12 lattice. 



self-energy. The resulting antinodal spectral functions 
for U = 5.2 are shown as dashed and dotted lines in Fig. 
Eta), respectively. The shaded region denotes the bounds 
in which one would expect the true spectrum, according 
to the opposite FS tendencies (with DQMC over- and 
DCA underestimating gaps at small cluster sizes) of both 
methods. Note that the remaining uncertainty is still 
significant and that the bounds are not rigorous, due to 
numerical noise and the difficulties of the MEM. 

Our unbiased estimate of Ax(u>), shown as solid line in 
Fig. Eta), reduces these uncertainties drastically: We find 
that the spectral weight at low frequencies (\oj\ < 0.3) is 
much smaller than predicted by DCA, but still signifi- 
cant (i.e., larger than the raw DQMC prediction). The 
true peak height at \u\ m 1.1 is close to the average of 
the DCA and DQMC predictions. At large frequencies 
| a; | > 1.5, we find excellent agreement with the earlier 
DQMC estimates^ which shows that the DQMC FS er- 
rors are small in this region (cf. Fig. [6]) and also veri- 
fies the procedures for analytic continuation; in contrast, 
DCA is still far off (at N = 64). 

Compared with the results for [7 = 4 presented in Fig. 
[6j our unbiased result [solid line in Fig.EJa)] shows much 
stronger PG characteristics, as is certainly expected at 
the stronger interaction U = 5.2. Spectra for a full range 
of temperatures at this interaction are shown in Fig.[8^b) 
for a 12 x 12 system; these results can directly be com- 
pared with the second column in Fig. [SJ Already at the 
highest temperature T = 0.50 [dotted line in Fig. ED>)], 
the spectral peak is much broader, i.e., more spectral 
weight has been shifted away from the origin than at 
U = 4. This tendency towards more insulating behavior 
remains at lower T: The peak-to-peak width is about 
twice as large as for [7 = 4. At T = 0.18 (dash-dotted 
line), no spectral weight can be resolved at |w| < 0.5, so 
that the PG looks numerically like a full gap. In addi- 
tion, the characteristic PG temperature is clearly shifted 
upwards, with a well-developed PG already at T = 0.28; 
the dependence of T* on U will be studied more broadly 
in subsection IIV CI 
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FIG. 9. (Color online) Unbiased local spectra A{ui) (first col- 
umn) and momentum-resolved spectra A^ioj) for U = 4 and 
k along the path through the Brillouin zone illustrated in 
Fig. []Ja) i the pseudogap opens with strong k dependence at 
T < T* « 0.20. A local maximum in the spectral density at 
oj 7^ (arrow) is indicative of spin-polaron physics. 

B. Evolution of pseudogap in full 
momentum-resolved spectral function 

So far, we have presented results which, for given pa- 
rameters U and T, are of a similar nature as those pre- 
viously discussed in the literature. The main advances 
of our study of nodal and antinodal spectra are (i) our 
elimination of the (enormous) finite-size bias inherent in 
raw results and (ii) our explicit analysis of temperature 
effects. We will now turn to fundamentally new results, 
namely spectra with full momentum resolution. 

Figure [9] shows unbiased momentum-resolved spectra 
Ak(w) throughout the whole Brillouin zone, along the 
path indicated in Fig. [TJa), at weak coupling {7 = 4 
and in a temperature range 0.18 < T < 0.28; in ad- 
dition, the left column contains the local spectra A(u>), 
corresponding to an average over all k. We have chosen 
a path rXM'MX that contains the irreducible portion 
XM' of the noninteracting Fermi surface (at half fill- 
ing). The inclusion of this subpath allows us to study 
the nodal-antinodal dichotomy continuously and in de- 
tail; more generally, all variations along this path (where 
eu = 0) arise from a k dependent self-energy, i.e., effects 
beyond DMFT. 

At T = 0.28 (first row in Fig. [9]), the local intensity 
maxima are unique at each k point and agree rather well 
with the noninteracting dispersion (dashed line), ex- 
cept for the edges oj > 4. A well-defined quasiparticle 
peak at uj w (especially sharp near k = M' and more 
diffuse at k = X) is consistent with a Fermi liquid de- 
scription. This picture changes at T = 0.20 (second row), 
when the spectrum splits at k sa X, i.e., a pseudogap 
opens at the antinodal point, while the rest of the spec- 
trum (at momenta with ^ 0) is essentially unchanged. 
The gap size decreases smoothly on the line X — > M'. 
Only at T < 0.18 (third row) the QP is destroyed also at 



FIG. 10. (Color online) Test of MEM accuracy at U = 4, 
T = 0.28: The local spectral function A(cj) (dashed line), cal- 
culated from the local Green function G(r), agrees well with 
the average Ak.{iS) = jj^2 k Ak(ui) (solid line). Both curves 
reveal a dip at u < 0.8 which is absent at this temperature in 
Am'( u ) (short-dashed line) and Ax(lu) (dotted line). 

k = M'; a PG then extends over all momenta. 

Compared to the strong temperature dependence along 
the path XM', the spectra appear nearly unchanged in 
the rest of the BZ. In particular, a sharp dispersive quasi- 
particle like band, indicated by an arrow in the top panel, 
evolves from the X point about half way towards the T 
point (and, equivalently by particle-hole symmetry, from 
the X point towards the M point). We interpret this 
feature, which is not accessible in conventional DQMC 
studies at FS, as the formation of a spin polaron band 
(arrow in Fig. [9|), with an energy offset at lower T indi- 
cating the magnetic exchange scale. It ends (at higher 
|ek|) in a "waterfall" which breaks up the band structure 
into low and high energy features.— 

Taken together, our results indicate that, apart from 
incoherent features at |w| ~ 4 which are present at all 
temperatures and should continuously evolve into Hub- 
bard bands with increasing U, interaction effects come 
into play with lowering T first very locally (in momen- 
tum space) around the antinodal X point; apparently, 
the strong enhancement of scattering by the van Hove 
singularity at X completely determines the physics in 
this region. This explains why the spectra can become 
sharper, implying a reduction in the imaginary part of 
the self energy, on the path from X towards T (up to 
the position of the arrow in Fig. [9l corresponding to the 
energy u indicated by dotted lines), i.e., with increasing 
£k and uj; a behavior which is exactly opposite to usual 
Landau Fermi liquid and also to DMFT physics. 

This suppression of spectral weight around X already 
at elevated temperatures also explains the slight dip seen 
in the local spectrum at T = 0.28 (dots and dotted lines 
in top panel of Fig. 02 cf. also Fig. IT0|) : While the mo- 
menta around M' and in the spin-polaron band region 
(arrow) contribute "normally" to the local spectrum, the 
contributions from momenta near X are spread out to 
about a much larger width (with a significant fraction at 
| a; | > 1); the missing weight at w <C 1 results in the dip. 

One might worry that this analysis puts too much con- 
fidence in the accuracy of our data and that the dip in 
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FIG. 11. (Color online) Spectral functions Ak(ui) for k along 
the Fermi edge (line from X — > M'; cf. Fig. [TJ. FS results 
(L — 8, 16) converge (slowly) to the thermodynamic limit 
both point-wise and by refinements of the k resolution. 

the local spectrum at T = 0.28, a local suppression in 
A(u) by abound 10% in a narrow frequency range, corre- 
sponding to a "missing weight" of about 1%, could also 
result from uncertainties in the MEM procedure. There- 
fore, we have checked its consistency and accuracy in the 
largest finite-size system (16 x 16) by comparing the lo- 
cal spectrum A(u) (dashed line in Fig. [10]), obtained by 
direct analytic continuation using MEM from the local 
Green function G(r) with the average of all (here 256) 
momentum- resolved spectra Ak(uj) in Fig. [TO] AsG(r) = 
-k J2k Gk( T ) = Gk(r), both spectra should agree, if eval- 
uated exactly: A(u) = Ak(u) = jj A^uj) . As the 
MEM is inherently nonlinear, due to the entropy con- 
straint, deviations must be expected in practice. How- 
ever, our procedure, with very accurate DQMC data, 
seems to be quite stable: Although the k dependent spec- 
tral functions A^u) differ substantially at different k 
points (shown in Fig. [10] only for the nodal and antinodal 
points using short-dashed and dotted lines, respectively) 
and have much more pronounced features than the local 
spectral function A(uS) (long-dashed line), their average 
Ak(uj) (solid line) agrees with it nearly within linewidth; 
only the magnified inset reveals tiny differences at small 
frequencies. So we conclude that our techniques are more 
than adequate and that the small dip discussed above is, 
indeed, physical. 

Let us, finally, stress that our eliminations of finite-size 
errors have been absolutely essential for obtaining unbi- 
ased momentum-resolved spectra, as illustrated in Fig. 
QT] for the path X — > M': Not only is the convergence at 
the end points k = X and k = M' slow, the k resolution 
is also quite coarse, with only one intermediate point for 
L = 8 and only three intermediate points for L = 16. 
It is clear that a very significant extension of the cluster 
size (e.g. to 64 x 64, implying a factor of 4 6 = 4096 in 
computer time) would be needed in order to match the 
momentum resolution of our extrapolation procedure. 

C. Evolution of characteristic pseudogap 
temperature T* with interaction U 

Apart from yielding a momentum dependent T* , the 
criterion used in subsection IIV Al has the disadvantage of 
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FIG. 12. (Color online) Properties of finite clusters and 
in the thermodynamic limit at U = 4: (a) Difference be- 
tween nodal and antinodal Green functions versus tempera- 
ture; its maximum defines the characteristic PG temperature 
T* ~ 0.20. (b) Unnormalized spin structure factor at anti- 
ferromagnetic wave vector k = (tt, it). Dotted vertical lines 
mark T* . 

depending on the ill-conditioned analytic continuation of 
the imaginary-time DQMC Green functions to the real 
axis. On the other hand, it is difficult to define spe- 
cific PG criteria on the level of the imaginary-time Green 
functions [cf. Fig. BJa) and BJb)] *2£ However, the nodal- 
antinodal dichotomy, i.e., the momentum dependence of 
the Green functions along the line X — > M' (arising from 
a momentum dependence of the irreducible self-energy) 
turns out to be illuminating: Fig. [T27 a) shows that the 
norm of the difference between the imaginary-time Green 
functions, 

\Gm> ~Gx\ = { f dr\G M ,{r) - G x (t)\ 2 / '/?} 1/2 , 

is strongly enhanced (at U = 4) in the temperature range 
where the PG opens. Not surprisingly, this peak be- 
comes sharper and shifts towards lower T in the ther- 
modynamic limit; the position of the maximum yields a 
natural unique definition of the characteristic PG tem- 
perature T* m 0.20, indicated by a vertical dotted line in 

Fig.[H 

As discussed in the introduction, the PG is associated 
(at n ks 1) with AF correlations and may be regarded 
as a precursor of a fully gapped long-range ordered AF 
phase which, in d = 2, is realized only in the ground 
state.— Thus, we should expect to see a strong enhance- 
ment in suitable spin correlation functions. While the 
nearest-neighbor spin correlations are only very moder- 
ately enhanced at T < T* (not shown), the spin structure 
function is seen in Fig. [L27 b) to increase by a full factor 
of 4 in the range 0.9 T* <T < LIT*. At the same time, 
FS effects explode at T < T*. All this shows that the 
PG is driven by the development of AF order at a scale 
which is large compared to the lattice spacing. 

The PG physics and, in particular, the momentum de- 
pendence observed at U = 4 should disappear at strong 
coupling, when already the high-temperature phase is 
gapped at n = 1^ The dichotomy should also vanish in 
the limit U — > 0, where the energy scale T spin vanishes, 
and so does the magnitude of the pseudogap. Indeed, 
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FIG. 13. (Color online) Difference between nodal and anti- 
nodal Green functions versus temperature: DQMC results at 
weak to intermediate coupling (2 < U < 8) for 12 x 12 clusters 
(thick lines) and 8x8 clusters (thin lines). Inset: associated 
T* in comparison with DMFT Neel temperature. 

the momentum dependence is seen in Fig. [13] to peak at 
U ~ 4 and to decay quickly for larger couplings, where 
also FS effects (which can be estimated from the thin 
lines, corresponding to 8 x 8, in comparison to the main 
12 x 12 results) become irrelevant. At fixed cluster size, 
also the results at weaker coupling (U = 3, U = 2) fall 
off; unfortunately, they suffer from significant FS effects 
which are too costly to eliminate. Still, the peak positions 
allow us to estimate T*(U) in the full range of weak to 
intermediate coupling as denoted by symbols in the inset 
of Fig. Q2J 26 Also shown is the mean-field estimate of the 
critical temperature for AF long-range order (solid line). 
At first sight, this DMFT estimate of the Neel tempera- 
ture T° MFT would appear irrelevant, as the true T N = 
by the Mermin- Wagner theorem. However, we find that 
T* re 0.9T£ MPT for 4 < U < 8; a correction of FS ef- 
fects for U — 2 and U = 3 should push the corresponding 
values of T* also below T^ MFT . So the DMFT identi- 
fies the relevant temperature scale for spin coherence (as 
was previously observed in the strong-coupling regime 2 ^) ; 



however, it lacks the momentum resolution which is es- 
sential to capture the pseudogap physics explored in this 
paper. 



V. CONCLUSION 

After decades of research, our understanding of the 
two-dimensional Hubbard model, especially regarding 
the extent to which it captures the pseudogap and high- 
T c physics of cuprates, is still far from complete. Numer- 
ical simulation a 21 ! 28 " — give valuable hints, but continue 
to be dominated by finite-size effects^ We have over- 
come the finite-size barrier and presented momentum- 
resolved spectral functions in the thermodynamic limit, 
obtained by systematic extrapolation of DQMC Green 
functions (L —> oo and At —> 0). Based on this achieve- 
ment, we were able to disentangle the delicate interplay 
of dynamical and spatial magnetic correlations. At weak 
to intermediate couplings, this interplay leads, indeed, to 
the formation of a pseudogap in the half-filled band. The 
pseudogap originates from a strong k dependence of the 
self-energy, which results in a d-wave-like anisotropy in 
the opening of the charge gap and a "waterfall" substruc- 
ture of the spectrum. The associated temperature scale 
T* is determined by the onset of antiferromagnctic fluc- 
tuations (and nearly agrees with T^ MFT ), i.e., is rather 
high compared to other coherence scales and should be in 
reach of experiments with ultracold fermions on optical 
lattices^^ 
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